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This paper outlines the methodology of interpreting channel 1 and 2 AVHRR radiance data over 
the oceans and describes a detailed analysis of the sensitivity of monthly averages of retrieved 
aerosol parameters to the assumptions made in different retrieval algorithms. The analysis is based 
on using real AVHRR data and exploiting accurate numerical techniques for computing single and 
multiple scattering and spectral absorption of light in the vertically inhomogeneous atmosphere- 
ocean system. We show that two-channel algorithms can be expected to provide significantly 
more accurate and less biased retrievals of the aerosol optical thickness than one-channel 
algorithms and that imperfect cloud screening and calibration uncertainties are by far the largest 
sources of errors in the retrieved aerosol parameters. Both underestimating and overestimating 
aerosol absorption as well as the potentially strong variability of the real part of the aerosol 
refractive index may lead to regional and/or seasonal biases in optical thickness retrievals. The 
Angstrom exponent appears to be the most invariant aerosol size characteristic and should be 
retrieved along with optical thickness as the second aerosol parameter. 
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1. Introduction 


The effect of tropospheric aerosols on global climate via the direct and indirect radiative forcings is 
one of the largest remaining uncertainties in climate change studies. 1 Current assessments of the 
direct aerosol radiative effect mainly focus on sulfate aerosols (e.g., Roeckner et al 2 ). It has 
become clear, however, that other aerosol types like soil dust and smoke from biomass burning and 
sea salt are also likely to be important climate forcing factors. 3 ’ 4 The magnitude and even the sign 
of the climate forcing caused by these aerosol types is still unknown. Once the globid distribution 
of aerosol properties such as the optical thickness, size distribution, and chemical composition is 
available, the calculation of the direct aerosol forcing using general circulation models (GCMs) is 
rather straightforward. 5,6 However, estimates of the indirect aerosol effect require information on 
the distribution of the aerosol number density and additional knowledge of the physics and 
chemistry of aerosol-cloud interactions, 7 which are still poorly understood. 

The retrieval of the global distribution of aerosol properties and determination of trends in its 
temporal variation can only be achieved using long-term satellite measurements. The standard 
one-channel Advanced Very High Resolution Radiometer (AVHRR) aerosol retrieval algorithm 8 
uses channel 1 radiances (nominal wavelength Xj = 0.65pm) and relies on the fact that the 
radiance reflected by an aerosol layer over the dark ocean surface is nearly proportional to the 
product of the aerosol phase function P(Q ) at the observation scattering angle 0, single- 
scattering albedo xn , and optical thickness x . The phase function and the single-scattering albedo 
are, in turn, dependent on the aerosol composition and size and shape distributions and are 
functions of many parameters. Even in the simplest case of a monomodal polydispersion of 
homogeneous spherical particles, the number of unknown model parameters is at least 5: aerosol 
optical thickness, real and imaginary parts of the refractive index, and effective radius and effective 
variance of the size distribution. 9 Since the AVHRR algorithm utilizes only one datum per pixel 
(channel 1 reflectance at a single observation geometry), it can retrieve only one model parameter 
(optical thickness), whereas all remaining parameters must be fixed a priori. Although the choice 
of the latter parameters can be optimized such that the algorithm produces a minimal long term 
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statistical deviation from existing sun photometer measurements of the optical thickness, 8 the 
strong spatial and temporal variability of the aerosol size, shape, and composition makes inevitable 
large errors in the assumed phase function and thus in the retrieved optical thickness in particular 
cases. Such large errors have been demonstrated not only by detailed sensitivity analyses, 10-12 but 
also by direct comparisons of AVHRR retrievals with sun photometer measurements. For 
example, figure 4 of Stowe et al . 8 shows relative errors in the retrieved x exceeding 100% and 
absolute errors exceeding 0.15. The results of Ignatov et a/. 13 show even larger discrepancies and 
suggest that more detailed comparisons of single-channel AVHRR retrievals with future sun 
photometer measurements are likely to reveal more significant errors. 

Another important limitation of the standard AVHRR algorithm is that it provides no 
information about the effective particle size and thus makes impossible estimates of the aerosol 
indirect radiative forcing. Indeed, the quantification of the Twomey effect 14 requires accurate 
satellite measurements of cloud condensation nuclei (CCN) column densities. The only way of 
retrieving the number of tropospheric aerosols in the vertical column of unit horizontal cross 
section from satellite radiance measurements is to divide the satellite-retrieved aerosol optical 
thickness x by the average extinction cross section per particle. Since the AVHRR algorithm 
assumes rather than retrieves the aerosol model, the strong sensitivity of the extinction cross 
section to assumed aerosol effective radius and the significant temporal and spatial aerosol 
variability make single-channel AVHRR retrievals of the CCN column concentration highly 
inaccurate. 15 

It has been suggested that the use of multichannel reflectance measurements can provide 
additional information on the aerosol model and also improve the accuracy of the optical thickness 
retrieval. 16-22 In addition to channel 1 radiance data, AVHRR provides channel 2 reflectances 
(nominal wavelength X 2 = 0.85pm ) that can be used to improve the performance of the AVHRR 
algorithm by retrieving two aerosol parameters rather than just one. Again, one of these 
parameters must be the aerosol optical thickness at a visible wavelength, whereas there is, in 
general, a choice for the second retrieved parameter. Indeed, even if the aerosol size distribution is 
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monomodal and the refractive index is wavelength-independent, one has a choice of retrieving the 
effective radius, the real part of the refractive index, the imaginary part of the refractive index, or 
the effective variance. The situation becomes even more complicated if the size distribution is 
bimodal or multimodal and/or if the refractive index varies with wavelength. 

For example, Nakajima and Higurashi 19 used a modified power law size distribution and 
retrieved as the second aerosol parameter the power exponent assuming that the aerosol refractive 
index is fixed. Higurashi and Nakajima 22 employed a bimodal log normal volume distribution, 
assumed the same fixed refractive index for both modes, and retrieved the relative contribution of 
modes 1 and 2 to the total aerosol number density. Obviously, other algorithms are possible. 16 It is 
thus clear that the performance of several candidate two-channel algorithms must be examined 
before a “standard” algorithm is selected as the one providing the best statistical accuracy. 

Another major issue is cloud screening. Even a small cloud contamination of a pixel, if not 
detected, can lead to a gross overestimation of the retrieved aerosol optical thickness. Several 
algorithms using different AVHRR spectral channels have been proposed (e.g.. Refs. 22-25). 
However, only the International Satellite Cloud Climatology Project (ISCCP) algorithm has been 
thoroughly validated, 26 and even this algorithm may need to be modified because it was primarily 
designed as a conservative cloud detection algorithm (most pixels for which the presence of a 
cloud is in doubt are declared cloud-free), whereas aerosol retrievals may need a more conservative 
cloud screening algorithm (pixels for which the presence of a cloud is in doubt are declared 
cloudy). 

One of the main objectives of the Global Aerosol Climatology Project (GACP; URL: 
http://gacp.giss.nasa.gov), established in 1998 as a joint initiative of NASA’s Radiation Science 
Program and Global Energy and Water Cycle Experiment (GEWEX), is to retroactively infer the 
global distribution of aerosols, their properties, and their seasonal and interannual variations for the 
full period of available satellite data. This is planned to be accomplished primarily through a 
systematic application of multichannel aerosol retrieval algorithms to existing satellite data and 
advanced 3-dimensional aerosol chemistry/transport models (e.g., Tegen et al. 2 ). In this paper we 
outline the methodology of interpreting channel 1 and 2 AVHRR radiance data over the oceans 
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and describe a detailed analysis of the sensitivity of retrieved aerosol parameters to the 
assumptions made in the retrieval algorithms. 

Given the strong spatial and temporal variability of tropospheric aerosols and the obvious 
limitations of two-channel algorithms, it is unreasonable to expect a high accuracy for each 
instantaneous retrieval. On the other hand, one may expect that imperfections of an algorithm 
may be partially compensated for by averaging the retrieved aerosol parameters over a 
sufficiently long period of time. In other words, although a two-channel algorithm cannot be 
expected to provide accurate daily retrievals, it may still provide rather accurate monthly, 
decadal, and annual averages that can be quite useful in climate research. Therefore, the strategy 
that we adopted for our sensitivity analysis was to work with real AVHRR data rather than with 
computer-generated, synthetic data and to look at the effect of various a priori assumptions made 
in specific candidate algorithms on monthly averages of the retrieved aerosol parameters. 

Another feature of our approach is the use of accurate numerical techniques for 
computing single and multiple scattering and spectral absorption of light in the vertically 
inhomogeneous atmosphere-ocean system. Two-channel satellite aerosol retrievals represent a 
complex underdetermined problem and make unavoidable many a priori assumptions regarding 
the parameters of the atmosphere-ocean model. Furthermore, the potentially strong 
contamination of AVHRR channel 2 radiances by water vapor absorption requires a special 
treatment. Therefore, the use of accurate numerical methods enabled us to focus on analyzing the 
effect of inherent uncertainty in model parameters rather than on an examination of possible 
retrieval artifacts resulting from the use of approximate solution approaches. 

2. Satellite Data, Atmosphere-Ocean Model, and Radiative Transfer Code 

Our initial activity has focussed on applying several aerosol retrieval algorithms to AVHRR 
channel 1 and 2 radiance data over the oceans contained in the gridded ISCCP DX dataset. The 
main advantages of using the ISCCP data product 28 are that it is easily available and contains an 
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elaborate cloud detection algorithm that can be easily modified for the purposes of aerosol 
retrievals. 

As an initial approximation, we assume that aerosol particles are homogeneous spheres and 
compute their scattering and radiative properties using the standard Lorenz-Mie theory. 9 The 
following aerosol parameters serve as an input for the single-scattering Mie code: the type of the 
aerosol size distribution, size distribution parameters, and real and imaginary parts of the 
refractive index. A single value is specified for each parameter except for the parameter to be 
retrieved. For the latter, a grid of values covering the expected range of its variation is provided. 
The output consists of a file containing the Legendre expansion coefficients of the aerosol phase 
function, the extinction cross section, and the single scattering albedo for the set of wavelengths 
used by the multiple-scattering code. 

Theoretical channel 1 and 2 reflectances are calculated using a multiple-scattering code 
based on the scalar version of the adding/doubling method. 9 The code takes into account the rough 
ocean surface reflection via the modified Kirchhoff approximation, 11 water vapor, oxygen, and 
CO 2 absorption via the k-distribution technique, and multiple scattering by stratospheric and 
tropospheric aerosols and molecules. 

The distribution of ocean surface slopes is assumed to be Gaussian, 

dz dz ' | _ 1 

dx’ dy) 2ns 1 

where the mean square surface slope s 2 is related to the near-surface wind speed W (m/s) via the 
empirical formula 30 

2s 2 =0.003 + 0.005 \2W. (2) 

The respective FORTRAN code 11 computes the Fourier components of the ocean bidirectional 
reflection function using the same grid of Gauss quadrature nodes representing the cosines of the 
angles of incidence and reflection as the adding/doubling routine. The number of Fourier 
components must be greater than or equal to the number of Fourier components used in the 
adding/doubling calculations for the atmosphere. Reflection geometries within 40° off the sun- 
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glint direction are excluded from the analysis. The upwelling radiances from the ocean body and 
foam scattering are either ignored or modeled by adding a small Lambertian component to the 
surface bidirectional reflection function. 

The gases that have lines in the first and second AVHRR channels are H2O, CO2, and O2, 
of which only H2O and O2 are significant contributors, although CO2 is included for the sake of 
completeness. The gaseous continua that are included are the ozone Chappuis band and the water 
vapor continuum . 31 The atmospheric temperature and humidity profiles are taken from the ISCCP 
version of the Television and Infrared Observation Satellite (TIROS) Operational Vertical Sounder 
(TOVS) data. The total atmosphere is subdivided into a number of homogeneous layers which is 
increased until convergent results are obtained. In most cases a ten-layer model provides quite 
sufficient accuracy. The vertical distribution of ozone and water vapor is based on a standard 
atmospheric profile . 32 Variations in ozone and water vapor are dealt with by scaling the total 
column amounts appropriately while maintaining the same normalized profile. The vertical 
profile of the aerosol number density is taken to be the same as the normalized profile of water 
vapor. Differences between this assumed profile and the actual profiles of water vapor and 
aerosol are not generally a significant error source in the radiative transfer modeling. There will 
be some errors in the retrieved aerosol size and optical depth when the majority of the aerosol is 
above the majority of the water vapor, or when the majority of the aerosol is below the majority 
of the water vapor. The chosen profile of aerosols provides a balance between these two 
extremes. Stratospheric aerosols are treated separately using aerosol size, optical depth and 
vertical profile information from Stratospheric Aerosol and Gas Experiment (SAGE) III . 33 

The radiative transfer code described above can be used to compute a look-up table for 
any candidate aerosol retrieval algorithm. Each look-up table is a file in which multidimensional 
arrays of theoretical channel 1 and 2 reflectance values for all viewing geometries and aerosol 
and atmospheric parameters are stored. The overall dimension of the table is determined by the 
product of the following parameters: 

• relative satellite-sun azimuth angle grid size; 
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• viewing zenith angle grid size; 

• solar zenith angle grid size; 

• ozone amount grid size (applies only to the channel 1 part of the look-up table); 

• water vapor amount grid size (applies only to the channel 2 part of the look-up table); 

• aerosol optical thickness grid size; 

• second retrieved aerosol parameter grid size. 

The look-up tables are used to retrieve the aerosol optical thickness and a second aerosol 
parameter using cloud-screened channel 1 and 2 radiance data. The retrieval routine performs a 
two dimensional search of the minimum of the respective error function using the so-called 
direction set method, 34 which does not require the calculation of derivatives. The error function 
is defined as 



where the subscripts 1 and 2 label scaled radiances in the first and second AVHRR channels and 
the letters t and m label theoretical and measured quantities, respectively. Iterations continue 
until the error function becomes smaller than a certain threshold value. Each pixel is mapped on 
a rxl* global grid. The retrieved values for all pixels within one grid cell are averaged to 
produce a map for a specified period of time. 

It is well known that multidimensional minimization may be a very complicated process and 
often may result in finding a local rather than a global minimum. Therefore, we have analyzed 
many particular cases “by hand” and made sure than in all cases considered the minimum found 
was the global minimum within the specified range of variability of model parameters. 
Furthermore, in addition to the direction set method, we have also implemented the so-called 
downhill simplex method. 34 The excellent agreement between the numbers obtained with the two 
quite independent minimization procedures may also indicate that our retrieval scheme produces 
reliable results. 
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3. Sensitivity Analysis 
A. Benchmark Atmosphere-Ocean Model 

Because the number of candidate algorithms is, in principle, unlimited, we decided to make the 
scope of our sensitivity analysis manageable by adopting the simple approach of selecting a 
benchmark atmosphere-ocean model and then examining the changes in the retrieved aerosol 
parameters caused oy variations in adopted model parameters. As the benchmark atmosphere- 
ocean model, we; selected the one based on a modified power law size distribution of the form 


n(r)- 
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with r x = 0.1pm, r 2 = 10pm, and a e[2.5, 5], The constant C is uniquely determined from the 
standard normalization 

00 

J dr n(r) - 1 . ( 5 ) 

o 

Note that larger values of the power exponent a correspond to smaller aerosols and vice versa. 
The refractive index is assumed to be wavelength-independent and equal to 1 .5 + 0.005/. 

Figure 1 shows the respective phase function versus scattering angle and power exponent a . 
A prominent phase function feature is a significant deepening of the valley at side scattering angles 
with decreasing a caused by the increasing effect of absorption inside larger aerosols. Figure 2 
plots effective radius r tff versus a , where 9 


J drrnr 2 n(r) 

r eff = ^ • (6) 

J dmr 2 n(r) 

0 

Figures 3 and 4 show the Angstrom exponent A as a function of a and r eff , respectively, where 
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A feflW] 
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( 7 ) 


lx=x, 

Finally, Fig. 5 shows the power exponent dependence of the single scattering albedo rn . ns 
remains nearly constant for a a 4, but significantly decreases for smaller a due to increasing 
absorption inside larger aerosol particles. 

The contribution of the upwelling radiation from the ocean body and foam scattering is 

q e 

ignored and the wind speed is fixed at a globally uniform value of W= 7 m/s. 

B. Cloud Screening 

Panels (a) and (b) of Fig. 6 show monthly mean optical thickness < x > and optical-thickness- 
weighted power exponent <a> for July of 1986 retrieved using the standard ISCCP cloud 
detection scheme, 24 the ISCCP calibration of the first AVHRR channel, 28 and the pre-launch 
calibration of the second channel. The ISCCP cloud detection scheme includes the following five 
major steps: (1) applying a space contrast test to individual infrared (IR) images; (2) applying a 
time contrast test to three consecutive IR images at constant diurnal phase; (3) cumulation of 
space/time statistics using both IR and visible images; (4) construction of composite clear-sky 
visible radiances and IR temperatures once every 5 days at each diurnal phase and location; and (5) 
applying radiance and temperature thresholds by comparing measured radiances and IR 
temperatures to the respective composite values. 

The relatively large optical thickness values in Fig. 6(a) may indicate a significant residual 
cloud contamination of many pixels that were classified as clear sky and suggest that the standard 
ISCCP criteria for detecting clear sky pixels may need to be tightened. It is well known that the 
effect of aerosols and clouds on visible channel reflectances is similar, whereas AVHRR channel 5 
(X 5 =1 1.7pm) reflectances are not affected by aerosols because of their negligibly small optical 
thickness at infrared wavelengths. Therefore, it is likely that tightening the visible ISCCP 
threshold or imposing an additional visible radiance threshold rejecting pixels with channel 1 or 
channel 2 radiances exceeding a certain value may result in an adverse loss of pixels with 
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significant aerosol loads. On the other hand, tightening the IR threshold can be expected to mostly 
affect only the results of cloud detection. Panels (c) and (d) of Fig. 6 present the retrieval results 
obtained with a modified ISCCP cloud screening scheme which retains only pixels with IR 
temperatures warmer than the composite values. Contrasting Figs. 6(a) and 6(c) shows a 
significant overall decrease in < x > , but a relatively weak effect on the pixels with large values of 
the optical thickness. Another noticeable effect is a significant overall increase in the mean power 
exponent and, thus, a decrease in the average particle size. Both effects are consistent with the 
assumption that the more conservative cloud screening algorithm removes more cloud- 
contaminated pixels and the fact that typical cloud particles are larger than typical aerosols. The 
results of Refs. 36-38 may indicate that the types of clouds eliminated by this algorithm are small 
cumulus clouds and optically thin cirrus. 

Panels (a) and (b) of Fig. 7 show the results obtained with an even more conservative cloud 
screening scheme, which retains only pixels wanner than the respective composite temperatures by 
IK or more. It is obvious that the overall < x > is further decreased, whereas the cases of large 
aerosol loads are hardly affected. The effect on < a > is weaker and less obvious. Large < a > 
seem to further slightly increase and small < a > seem to slightly decrease, but this may be an 
artifact of reducing statistics (see below). 

Wagener et al 25 suggested to use the channel 1 to channel 2 radiance ratio, S i2 , as an 
additional indicator of cloud contamination by rejecting all pixels with S n <15 or S n >3.5. 
This criterion is based on the observations that S 12 £ 1 for totally overcast pixels or data over land 
surfaces. Panels (c) and (d) of Fig. 7 show the results obtained by superimposing the S n criterion 
on the more conservative IR threshold used for Figs. 7(a) and 7(b). The overall change is 
relatively small, although some spurious cases of large x at high northern and southern latitudes 
are removed. 

The results of the above sensitivity tests clearly demonstrate that accurate cloud screening is 
an issue of critical importance. Applying increasingly conservative thresholds may further reduce 
the risk of cloud contamination, but also reduces the amount of useful aerosol data [as already 
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manifested by the increased number of white pixels in Fig. 7(c) as compared to Fig. 6(a)] and may 
ultimately introduce a significant statistical bias. For example, the latter three cloud screening 
schemes rejected, respectively, 56.5%, 78.8%, and 79.9% of the pixels originally classified by 
ISCCP as cloud-free. The latter two numbers are consistent with the estimate of the amount of 
clouds missing by ISCCP derived by Liao et al 31 from SAGE II data. Therefore, we decided to 
adopt for the following analysis the combination of the conservative IR scheme (retaining only 
pixels warmer than the composite values by IK or more) and the 1.5<S'i 2 <3.5 criterion. It is 
obvious, however, that the definitive examination of the quality of the product generated by this 
algorithm will require extensive comparisons with long-term ground-based measurements and, 
possibly, future results from more advanced satellite instruments. 39 

The most obvious features of the optical thickness patterns in Fig. 7(c) are the plumes of 
African and Asian dust and aerosols produced by biomass burning in Equatorial Africa and South 
America. Dust aerosols can also be clearly identified in the power exponent map [Fig. 7(d)] due to 
their larger sizes. Figure 7(d) shows a remarkable asymmetry in the average aerosol size: the 
northern hemisphere seems to be dominated by relatively small particles, presumably 
anthropogenic pollutants, whereas a large fraction of the southern hemisphere is covered by 
significantly larger particles, most likely sea salt aerosols. Figure 8 depicts the global average as 
well as the northern and southern hemisphere averages of the aerosol optical thickness and optical- 
thickness-weighted power exponent retrieved over the full period of NOAA-9 observations. The 
aerosols in the southern hemisphere indeed appear to be systematically larger as well as optically 
thinner than those in the northern hemisphere. 

C. Effect of Radiance Calibration Uncertainties 

The retrievals described in the previous section were based on the ISCCP post-launch calibration 
of the NOAA-9 AVHRR channel 1 and the pre-launch calibration of channel 2. It is known that 
the in-flight degradation of channel 2 was significantly slower than that of channel 1 . 40 This may 
explain why neither curve in Fig. 8 shows a significant long-term trend. There may be a slight 
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overall increase in the power law exponent, which would be consistent with the fact that the in- 
flight degradation of channel 1 was and that of channel 2 was not corrected for, thereby 
increasing the spectral contrast between the channels and reducing the retrieved particle size. 
However, even if this trend is real, it is very weak. 

The post-launch calibration of AVHRR channels 1 and 2 is a very complicated problem 
associated with many uncertainties and discrepancies. 40,41 To examine the potential effect of 
calibration uncertainties on the two-channel aerosol retrievals, we have recomputed Figs. 7(a) 
and 7(b) using the NOAA post-launch calibration of both visible channels. 40 Comparison of Figs. 
7(a) and 7(b) with Figs. 9(a) and 9(b), respectively, shows significant changes in both the optical 
thickness and the power law exponent patterns, especially in the regions with small aerosol 
loads. We have found that these changes are mainly caused by very small differences in the 
assumed values of the so-called deep space count (the response of the radiometer to zero incident 
intensity) between the ISCCP and NOAA calibrations of channel 1 and the pre-launch and 
NOAA post-launch calibrations of channel 2. Although these differences cause percent scaled 
radiance differences as small as a few tenths of a percent, the effect on the accuracy of retrieving 
small optical thickness values is rather strong. Because aerosol retrievals over the ocean surface 
deal with small measured radiances and because the aerosol contribution to the total radiance is 
often weak, we have to conclude that the post-launch calibration issue for channels 1 and 2 may 
need to be revisited before one attempts a massive aerosol retrieval for the full period of 
available AVHRR data. 

There is no doubt that the two-channel AVHRR retrievals will be very useful for studying 
monthly and seasonal variability of spatial patterns of aerosol parameters. However, the very 
methodology that has to be used for the post-launch calibration 40,41 and its rather poor accuracy 
make it difficult to expect that these retrievals by themselves can detect a slow trend in the 
aerosol optical thickness and/or size with accuracy needed for climate change studies. 5 Our 
results demonstrate once again that reliable detection of a long-term trend in climatically 
important aerosol parameters may require the use of a much more accurate remote-sensing 
technique such as high-precision polarimetry. 11,12,15 
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D. Different Ways of Averaging the Power Exponent 

Figure 7(d) shows the optical-thickness-weighted monthly mean power law exponent computed as 
1 ,+T 

<a >=— f a(r)x(r) dr , (8) 

T <x > J 

t 

where 

1 ,+T 

<T>=— Jx(0d». (9) 

t 

This quantity is a measure of the average size of aerosols suspended in the atmosphere at a given 
location during the month. An alternative definition is 

1 ,+T 

a = -J a (0d' (10) 

r 

and indicates the aerosol size encountered most frequently during the same period of time without 
indicating how much aerosol had that size. < a > may differ significantly from a if there is a 
strong correlation between a and x , but should be nearly the same if the correlation is absent or 
weak. 

Figure 9(c) shows a map of a for July of 1986 and should be contrasted to Fig. 7(d), 
whereas Fig. 9(d) maps the ratio a / < a > . One sees that in most cases the difference between a 
and <a > is within ±10%, which indicates a weak correlation between a and x on the global 
scale. This conclusion is corroborated by Fig. 10, which presents a regression of a against x and 
shows no obvious correlation. On the other hand, Fig. 9(d) may indicate the presence of significant 
local correlations, especially in the southern hemisphere. 

E. Effect of Aerosol Absorption 

Recent studies have shown that a significant fraction of tropospheric aerosols (especially the 
mineral dust and biomass burning components) can be rather strongly absorbing. 42 ' 45 This may 
justify the choice of a nonzero imaginary part of the refractive index for a unified aerosol model 


14 



used in global AVHRR retrievals. However, many remote areas can be dominated by 
nonabsorbing aerosols such as sea salt, and the total single scattering albedo can be significantly 
closer to unity than is shown by the solid curve in Fig. 5. Panels (a) and (b) of Fig. 1 1 show the 
ratios of the average optical thickness and optical-thickness-weighted power exponent retrieved 
with the same aerosol model but assuming an imaginary part of the refractive index of Im(/w) = 
0.002 relative to those displayed in Figs. 7(c) and 7(d), respectively. The most obvious result of 
decreasing absorption is an overall decrease in the retrieved optical thickness. The decrease is 
especially significant and can exceed 25% in areas dominated by larger aerosols (smaller power 
exponents) as well as in areas with heavy aerosol loads. These changes can be explained by 
increased phase function values at side- and back-scattering geometries for larger particles (cf. 
Figs. 1 and 12) and by systematically larger single scattering albedos (cf. solid and dashed curves 
in Fig. 5). On the other hand, decreasing absorption does not seem to have a significant effect on 
the retrieved aerosol size [Fig. 1 1(b)]. 

It is thus clear that the use of a globally-unified imaginary part of the refractive index can 
result in significant systematic regional and/or seasonal errors in the retrieved aerosol optical 
thickness. Since Im(m) cannot be inferred from AVHRR channel 1 and 2 data, this result may call 
for applying several aerosol models with different Im(/w) values and using auxiliary information 
such as Total Ozone Mapping Spectrometer (TOMS) data 46 or results of aerosol transport 
modeling 3 as indicators of the likely magnitude of aerosol absorption. Of course, this would result 
in a much more complicated and time-consuming retrieval algorithm. Furthermore, absorption can 
be spectrally dependent, thereby potentially affecting the retrievals of the aerosol size. 

F. Effect of the Real Part of the Refractive Index 

It is well known that the real part of the aerosol refractive index may be highly variable in space 
and time and may differ significantly from the adopted benchmark value 1.5. 47 To examine 
potential retrieval errors caused by this variability, we performed retrievals similar to those 
shown in panels (c) and (d) of Fig. 7, but assuming a refractive index value of 1.4 + 0.005/. 
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Figures 13(a) and 13(b) show the ratios of the aerosol optical thickness and power exponent thus 
obtained relative to those shown in Figs. 7(c) and 7(d), respectively. The most notable result is a 
substantial overall increase of the optical thickness and a very little change in the power 
exponent. The former is obviously caused by systematically lower phase function and single 
scattering albedo values for Re(m) = 1.4 (cf. Figs. 1, 5, and 14). This test suggests that adopting 
a fixed, globally uniform refractive index may result in significant instantaneous retrieval errors 
and, potentially, in a systematic regional and seasonal bias in areas dominated by a single aerosol 
type with refractive index significantly different from the adopted value. A feasible, but 
cumbersome solution of this problem could be to constrain the refractive index range by 
employing auxiliary information provided, e.g., by aerosol transport models. 

G. Effect of Diffuse Ocean Reflectance 

Previous results were obtained assuming no upwelling radiance contribution caused by scattering 
beneath the ocean surface and by foam. Depending on meteorological conditions and location, this 
diffuse contribution can be significant 48-51 but is difficult to parameterize for the purpose of the 
global aerosol retrieval. Most of existing parameterizations require (often unavailable) information 
on the real-time wind speed and pigment concentration, were derived under specific natural or 
even artificial conditions, and may not be readily generalized. One may expect that the use of the 
previously described S 12 threshold may eliminate the cases of most significant foam 
contamination caused by wind speeds exceeding 12 m/s because foam tends to reduce the contrast 
between channel 1 and 2 radiances. Furthermore, by selecting relatively warmer pixels our 
algorithm favors the cases with higher ocean surface temperatures and, thus, lower wind speeds. In 
addition, pixels with very high wind speeds are more likely to be cloudy. For lower wind speeds, 
one may have to use the simplest parameterization of the diffuse ocean reflectance by adding a 
small, uniform Lambertian component to the ocean surface reflection function as suggested by 
Stowe et al* 
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Panels (c) and (d) of Fig. 1 1 demonstrate the effect of adding a fixed 0.003 Lambertian 
component that slightly exceeds the 0.002 component used in Ref. 8. The most obvious and 
natural result is an overall reduction of the optical thickness. This change is especially noticeable 
in the cases of low aerosol loads, but is within 10% in the majority of other cases. A secondary 
(and mostly negligible) effect is a slight increase in the retrieved power exponent caused by the 
implicit increase of the spectral contrast in the aerosol scattering contribution to the total 
reflectance. The latter effect may further weaken if the foam contribution itself shows a spectral 
contrast comparable to that of aerosols. 49,50 

H. Effect of Wind Speed 

The deviation of the actual wind speed from the benchmark 7 m/s value may also cause retrieval 
errors via changing the specular component of the ocean reflection function [Eqs. (1) and (2)]. 
Because the actual real-time value of the wind speed is often unknown and cannot be inferred from 
AVHRR data, one must examine the potential effect of adopting a globally uniform and constant 
wind speed value. Panels (c) and (d) of Fig. 13 computed for W = 11 m/s as well as similar 
computations for wind speed values smaller than 7 m/s (not shown) suggest that the retrieval errors 
in the optical thickness should be less that 10% in most cases, whereas the errors in the power 
exponent are negligibly small. Only in cases of very low aerosol loads ( x £ 0.05) the errors in the 
retrieved optical thickness can exceed 25%. These results are in a good agreement with those 
derived by Higurashi and Nakajima. Despite this relatively weak sensitivity, one may expect that 
making use of advanced global ocean wind speed data sets can further improve the accuracy of 
satellite aerosol retrievals. 

I. Effect of Size Distribution Function 

The modified power law used in the benchmark model is not necessarily the best representation of 
the actual shape of the aerosol size distribution in many cases, and other distribution functions such 
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as gamma and log normal distributions have often been used. Figure 1 5 summarizes the retrieval 
results obtained with a modified bimodal log normal distribution of the form 22 
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with r gl =0.17pm, r g2 =3.44jim, a gl =1.96, a gl =2.37, and ye[0. 1,100] and assuming the 

benchmark refractive index value m= 1.5 + 0.005/'. Again, the constant C] is determined from the 
normalization condition of Eq. (5). Figures 16 and 17 show the respective phase function and 
single scattering albedo. 

A very important conclusion following from Fig. 15(a) is that despite the large difference in 
the functional form of the power law and bimodal size distributions, the respective monthly mean 
optical thicknesses are remarkably similar and in most cases are within ±10% of each other. The 
retrieved monthly mean effective radii are in a much worse agreement and may differ by more than 
a factor of two. This may be explained by a significantly narrower range of possible r efr values for 


the bimodal size distribution and the saturation of A with increasing r efT for the power law 
distribution (Fig. 4). As a result, the retrieval algorithm based on the bimodal size distribution 
tends to produce smaller effective radii than that based on the power law distribution. The 
differences in the respective optical-thickness-weighted mean Angstrom exponents <A> as well as 
in the simple mean Angstrom exponents A [cf. Eqs. (8) and (10)] appear to be less significant 
than those in r eff . In most cases the differences in <A> and A do not exceed 0.3 and 0.2, 

respectively. These results may suggest that A is the aerosol size characteristic least sensitive to 
the uncertainties in the atmosphere-ocean model and should be retrieved along with optical 
thickness as the second aerosol parameter. 


J. Two-Channel Versus One-Channel Aerosol Retrievals 

In order to compare the expected performance of two-channel and one-channel retrieval 
algorithms, we repeated retrievals using the benchmark aerosol model with five fixed power 
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exponent values (a=2.5, 3, 3.5, 4, and 4.5), utilizing only AVHRR channel-1 radiance data, and 
retrieving only the aerosol optical thickness. The results of this test are summarized in Figs. 1 8 and 
19, which show a remarkable zonal structure most likely caused by seasonality of predominant 
AVHRR scattering geometries and systematic phase function differences for different values of a 
(Fig. 1). The differences in the retrieved optical thickness for different fixed a can exceed 300% 
and are indicative of biases that can be expected if a one-channel algorithm employs a fixed 
aerosol size not representative of the seasonally and/or regionally dominant aerosol type. As a 
whole, contrasting Fig. 15(a) with Figs. 18 and 19 convincingly demonstrates the great advantage 
of utilizing AVHRR channel 2 data in addition to channel 1 radiances. 

4. Conclusions 

The main results of our sensitivity study of the expected performance of two-channel aerosol 
retrieval algorithms based on real AVHRR data can be summarized as follows. 

• Two-channel algorithms can be expected to provide significantly more accurate and less biased 
retrievals of the aerosol optical thickness than one-channel algorithms. 

• Imperfect cloud screening and calibration uncertainties are by far the largest sources of errors in 
the retrieved optical thickness. Both problems are difficult to solve definitively and should be 
addressed by means of extensive ground-based observations, careful statistical analyses of the 
radiance data, and, potentially, comparisons with future results from more advanced satellite 
instruments. 

• Two different ways of computing the average aerosol size (direct versus optical-thickness- 
weighted) can be expected to produce similar results because of weak correlation between the 
aerosol optical thickness and size. 
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• Both underestimating and overestimating aerosol absorption as well as the potentially strong 
variability of the real part of the aerosol refractive index may lead to regional and/or seasonal 
biases in the retrieved aerosol optical thickness. 

• Deviations of the actual wind speed from the global mean value 7 m/s within the range from 0 
to 1 1 m/s have little effect on the retrieved optical thickness. Neglecting the diffuse component 
of the ocean reflection function can affect the retrieved optical thickness in the cases of low 
aerosol loads. 

• Simple monthly average of the Angstrom exponent appears to be the most invariant aerosol 
size characteristic and should be retrieved along with optical thickness as the second aerosol 
parameter. 

It may be expected that the two-channel retrieval algorithms could be improved by adopting 
time-dependent regional aerosol models. For example, dust particles have distinctly nonspherical 
shapes and their scattering properties can differ substantially from those for surface- or volume- 
equivalent spheres. 53-57 It is feasible, therefore, that using a nonspherical aerosol model can 
improve the accuracy of optical thickness retrievals over areas where dust aerosols are the 

CO 

dominant component. In addition, TOMS data could be used to determine areas with 
predominantly absorbing aerosols. 46 

Although our analysis demonstrates the expected range of retrieval errors caused by 
unavoidable uncertainties in the assumed parameters of the atmosphere-ocean model, it cannot 
determine the best choice of fixed, globally uniform values for all model parameters other than two 
parameters that are being retrieved. This choice can only be made on the basis of an extensive 
validation versus statistically representative ground-based and in situ measurements. This work 
has already begun (e.g., Refs. 22 and 59) and will be an important part of the GACP. Another 
aspect of sensitivity studies that needs to be addressed is the analysis of the effect of radiance 
calibration and water vapor amount errors on the retrieved aerosol parameters. It is expected that 
candidate retrieval algorithms will be further refined by using retrievals from POLDER 59 and 
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results from future space missions 39 and then re-applied to the full AVHRR data set. We also plan 
to use future airborne results from the Research Scanning Polarimeter built by SpecTIR 
Corporation . 60 

We thank two anonymous reviewers for useful comments on an earlier version of this paper 
and Nadia Zakharova for help with color graphics. This research was sponsored by the NASA 
Radiation Science Program managed by Robert Curran. 
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Figure captions 


Fig. 1. Phase function versus power exponent and scattering angle for the power law size 
distribution of Eq. (4) and m = 1.5 + 0.005/. 

Fig. 2. Effective radius versus power exponent for the power law size distribution of Eq. (4). 

Fig. 3. Angstrom exponent versus power exponent for the power law size distribution of Eq. (4). 

Fig. 4. Angstrom exponent versus effective radius for the power law size distribution of Eq. (4) 
and the bimodal size distribution of Eq. (11). 

Fig. 5. Single scattering albedo versus power exponent for the power law size distribution of Eq. 
(4) and three values of the refractive index. 

Fig. 6. (a) and (b): Monthly mean optical thickness < x > and optical-thickness-weighted power 
exponent <a > for July of 1986 derived using the benchmark atmosphere-ocean model and the 
standard ISCCP cloud detection scheme, (c) and (d): As in panels (a) and (b), but using a 
modified cloud detection scheme that retains only pixels with channel 5 temperatures warmer 
than the respective composite values. 

Fig. 7. (a) and (b): As in panels (a) and (b) of Fig. 6, but using the cloud detection scheme that 
retains only pixels with channel 5 temperatures warmer than their composite counterparts by 1 K 
or more, (c) and (d): As in panels (a) and (b), but with the addition of the 5,2 threshold. 

Fig. 8. Global and hemisphere averages of the aerosol optical thickness and power exponent 
plotted with daily resolution. The horizontal ticks correspond to the 15th day of the respective 
months. 

Fig. 9. (a) and (b): As in panels (a) and (b) of Fig. 7, but using the NOAA post-launch 
calibration of NOAA-9 AVHRR channels 1 and 2. (c): As in panel (d) of Fig. 7, but for a . (d): 
The a/ < a > ratio. 
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Fig. 10. Scatterplot of aerosol power exponent versus optical thickness for July of 1986. 

Fig. 11. (a) and (b): Ratios of the mean optical thickness and power exponent retrieved assuming 
the refractive index m =1.5 + 0.002/ relative to those shown in Figs. 7(c) and 7(d), respectively, 
(c) and (d): Ratios of the mean optical thickness and power exponent retrieved assuming a small 
added diffuse contribution to the ocean reflection function relative to those shown in Figs. 7(c) 
and 7(d), respectively. 

Fig. 12. Phase function versus power exponent and scattering angle for the power law size 
distribution of Eq. (4) and m - 1 .5 + 0.002/. 

Fig. 13. (a) and (b): Ratios of the mean optical thickness and power exponent retrieved assuming 
the refractive index m =1.4 + 0.005/ relative to those shown in Figs. 7(c) and 7(d), respectively, 
(c) and (d): Ratios of the mean optical thickness and power exponent retrieved assuming a wind 
speed of W= 11 m/s relative to those shown in Figs. 7(c) and 7(d), respectively. 

Fig. 14. Phase function versus power exponent and scattering angle for the power law size 
distribution of Eq. (4) and m- 1.4 + 0.005/. 

Fig. 15. (a): Ratio of monthly mean aerosol optical thicknesses retrieved using the bimodal size 
distribution of Eq. (11) and the modified power law of equation (4). (b): As in (a), but for the 
effective radius ratio, (c): As in (a), but for the difference of optical-thickness- weighted monthly 
mean Angstrom exponents, (d): As in (c), but for the difference of simple monthly mean 
Angstrom exponents. 

Fig. 16. Phase function versus y and scattering angle for the bimodal log normal size distribution 
of Eq. (11) and m = 1.5 + 0.005/. 

Fig. 17. Single scattering albedo versus y for the bimodal log normal size distribution of Eq. (11) 
and m = 1.5 + 0.005/. 
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Fig. 18. Aerosol optical thickness retrieved with the one-channel algorithm based on the power 
law size distribution of Eq. (4) with a= 2.5, 3, 4, and 4.5 relative to that retrieved for a= 3.5 for 
July of 1986. 

Fig. 19. As in Fig. 18, but for January of 1986. 
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